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Critical exponents of the infinitely slowly driven Zhang 
model of self-organized criticality are computed for d — 2,3 
with particular emphasis devoted to the various roughening 
exponents. Besides confirming recent estimates of some ex- 
ponents, new quantities are monitored and their critical ex- 
ponents computed. Among other results, it is shown that the 
three dimensional exponents do not coincide with the Bak, 
Tang, and Wiesenfeld (abelian) model and that the dynam- 
ical exponent as computed from the correlation length and 
from the roughness of the energy profile do not necessarily 
coincide as it is usually implicitly assumed. An explanation 
for this is provided. The possibility of comparing these results 
with those obtained from Renormalization Group arguments 
is also briefly addressed. 



I. INTRODUCTION 

Despite more of a decade of intensive studies, the phe- 
nomenon named Self-Organized Critically (SOC) by Bak, 
Tang, and Wiesenfeld (BTW) 0] is far from being fully 
understood. The name SOC originates from the funda- 
mental property that an open system externally driven 
in a (infinitely) slow fashion, settles into a critical state 
with no characteristic time and length scales, without 
any parameter tuning; see e.g. Rcf. pj for a review. 

Although a large number of recipes have been pro- 
posed as toy models to mimic this behavior, the origi- 
nal sandpile model [Q still carries most of the informa- 
tion presented in this phenomenon. A variation of this 
model was introduced a couple of years later by Zhang 
H . The basic differences with respect to the BTW model 
were: first, the variable describing the state of the lattice 
site could take continuous rather than discrete values; 
second, the BTW model is abelian Q while the Zhang 
model is not. In spite of this differences, extensive re- 
cent numerical simulations 0] on the two-dimensional 
Zhang model, opened the possibility that they both be- 
long to the same universality class, in disagreement with 
the original scaling prediction by Zhang |3|. Apart from 
the aforementioned investigation [^J, the Zhang model 
was already studied in different dimensionalities in Ref. 
P| where an estimates for some critical exponents, no- 
tably the avalanche size exponent r s , was given. How- 
ever, these estimates, whose main aim was to test the 
robustness of universality of the model under anisotropy 
of the energy repartition, appeared to be based on small 



sizes and statistics. 

On the other hand, a Langcvin counterpart of the 
Zhang model was repeatedly studied by Renormalization 
Group (RG) methods (?|,||,|9|,|l0| , and predictions for criti- 
cal exponents in a one-loop working scheme were drawn. 
The dynamical exponent z as calculated from the corre- 
lation function in the case when the additive noise has a 
typical time scale much bigger than the relaxation time 
scale ||, turned out to be very close to the one relat- 
ing the correlation length and the relaxation time in the 
standard dynamical scaling hypothesis [[tl] in the Zhang 
model. 

It is then desirable to have a more complete numer- 
ical investigation touching upon those issues appearing 
in the RG calculations and those which were previously 
neglected. This is indeed the aim of the present work 
where a fairly complete analysis of the model, in differ- 
ent dimensionalities is carried out and compared, when 
possible, with previous numerical and RG work. By do- 
ing this we found few unexpected results. 

Firstly, the three-dimensional results do not support 
the conjecture that the Zhang and the BTW models be- 
long to the same universality class. Secondly, whereas it 
is true that the exponent z of the Zhang model is very 
close to the one obtained by RG techniques as previously 
discussed, the roughening exponent is not ||. Finally, 
the critical exponent z is different when calculated from 
the dynamical scaling ansatz and when computed from 
the roughness exponent. This latter discrepancy can be 
fixed in our case by noting that the correlation length 
(maximum avalanche distance) does not scale linearly 
with system size L. 

The plan of the paper is as follows. In Sec. || the 
model is defined, whereas in Sec. HI all relevant quan- 



tities concurring to identify the critical behavior of the 
model are laid down. Sec. IV contains the results of this 



effort and comparisons with earlier ones, and finally some 
conclusive remarks are left in Sec. [V[ 



II. THE SLOWLY DRIVEN ZHANG MODEL 

Each point of an hypercubic lattice is characterized 
by a continuous energy variable E T (x, t) , where x de- 
notes the lattice position, t the driving (slow) time, and 
r the relaxation (fast) time. Whereas t runs from to 
a sufficiently large value needed to get good statistics, 
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r runs from to T(t), which is the total fast time that 
an avalanche initiated at a slow time t takes to be com- 
pleted. In this way the two time scales are well separated. 
Starting from an initially empty lattice, the dynamics of 
the evolution is defined as follows ||: 

1. Start with a randomly chosen lattice point Xo and 
set it slightly above some critical energy E c (here- 
after chosen to be 1 without loss of generality) by 
repeated addition of a random energy taken uni- 
formly from the interval (0, 1/4) fj^Jl3]] . 

2. The site xo relaxes according to the following equa- 
tion: 

E T+1 (x,t) = [1 - 0(Er(x,t) - E c )]E T (*,t) + (1) 

y(x) 

where #(•) is the Heaviside step function and d is the 
space dimension. Here the notation S y ( x ) means 
that the sum is restricted to the nearest-neighbors 
y of site x. Clearly this is tantamout to say that 
each site x whose energy exceeds a critical value E c 
is set to zero and its energy is equally redistributed 
to the nearest- neighbors. 

3. Iterate Step 2 for the other sites that become crit- 
ical until all sites are below E c . 

4. At this point increase t by one unit (t — > t + 1), and 
randomly pick a new initial seed Xq in Step 1. 

The process is iterated until the system has reached a 
steady-state configuration where the average energy 

X 

reaches a well defined value. Here V = L d is the volume 
of the lattice. We note that whenever there is no sub- 
script for the energy it will be implicitly assumed that 
the avalanche is over, i.e. r has reached T(t). Starting 
at this time, when the system has reached a stationary 
state, we collect all the relevant dynamical properties. 

III. PROBABILITY DISTRIBUTIONS AND 
CORRELATION FUNCTIONS 



T(t) 

S(t) = X>(t). (4) 

r=l 

From the size of the avalanche we can compute a charac- 
teristic length £(£) defined as the radius of gyration with 
respect to the seed site xo. This characteristic length is 
related to the time the avalanche needs to be completed 
through the standard relation 

T(t) (5) 

which defines the dynamical exponent z. 

Other quantities that are interesting to measure are 
the total input and output currents. They are defined as 
follows 

J in (i) =6E(xo,t) (6) 

T(t) 

r=OxGOA 

where A is the bulk and dA is the boundary of the bulk 
(the sum of the two forming the total available lattice 
space). Here SE(xo,t) is the total added energy neces- 
sary to the site to be active (i.e. above the critical energy 
E c = 1). 

In order to take into account the existence of two differ- 
ent time scales one should be very careful when defining 
the correlation functions. Upon extending (^), we can 
define the q— th spatial moment of the energy as 



£'W = 7 E^ f ) ( g ) 

X 

and then the interface width (or roughness) is 

W s {t,L) = y/W(f)-E(ff. (9) 



This definition applies to the slow time scale as also indi- 
cated by the suffix s, and coincides with the usual defini- 
tion of roughness in the framework of growth processes. 
On the other hand one could think to measure the energy 
fluctuations during the evolution of an avalanche. Since 
an avalanche of duration t occurs at many different input 
times t, we define the following fast roughness: 

W](t,L) = 



At each time t there is a growing avalanche; within the 
fast time scale we can measure the number of active sites 
at each update (t) 

S T (t) = J2 e ( E A*,t)-E c ) (3) 

X 

and from this we can define the size of an avalanche at 
time t 



In Eq. ([To]) the roughness is averaged over different times 
t (and hence avalanches). 

According to standard scaling hypothesis (see e.g. [fl4"|]) 
one expects these correlation functions to display the fol- 
lowing scaling forms: 
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Wf(r,L) = T - p t&f{r/L'f) 



W s (t,L)=t^<£ s (t/L z °) 



where anf $ s are finite size functions. 



(11a) 



(lib) 



In (11a) the roughness is expected to decrease rather 
than to increase as in more conventional growth processes 
PI , because the maximum energy is bounded and the 
avalanche is a relaxational process. 



IV. CRITICAL EXPONENTS AND RESULTS 

This model was already carefully investigated in two- 
dimensions. Apart from the original work 0], recent ex- 
tensive simulations on remarkably large sizes were carried 
out in d = 2 ||. When comparable, our results are in 
good agreement with both previous analysis. However, 
in these papers, the behavior of some important quanti- 
ties, necessary to our purposes, were neglected, nor was 
a complete study in dimensionality different from 2 ever 
attempted. Indeed while Zhang only reported the steady 
state value of the average energy along with the "quan- 
tized" energy distribution P(E) for d = 3, in fi a value 
for the avalanche size exponent (see below) is reported 
for dimensions up to 4. The latter was however probably 
based on very small sizes without any attempt of a finite 
scale analysis. As a result these estimates, albeit close, 
turn out to be slightly off compared to ours. 

In our simulations we used sizes up to L = 300, 60 and 
times up to 2 17 ,2 18 in d = 2,3 respectively. These are 
smaller than the ones used in || for d = 2 but consider- 
ably larger than all other three dimensional studies. 

For the sake of clarity and compactness, let us now 
review some known results first. As it is well known by 
now, the system reaches a steady state (where the aver- 
age energy is no longer changing) after a transient which 
clearly scales as L d since it takes that many time steps 
(on average) to "explore" the whole lattice. The result- 
ing values of the stored energy E(t) are 0.63 ± 0.01 and 
0.58 ± 0.01 (estimated from the largest sizes) in d = 2, 3 
respectively. These results are in agreement with those 
found by Zhang in his original simulations. Another fea- 
ture already observed by Zhang is that the critical state 
has energy which is peaked around well defined energies 
the number of which depends only on the dimensional- 
ity of the hypercubic lattice. It has also been established 
that this feature is unchanged upon introducing an asym- 
metry in the probability distribution and by introducing 
different lattices ||. 

As explained by Hwa and Kardar |l5| in the framework 
of the one-dimensional BTW sandpile model, monitoring 
the total output energy current proves to be very useful 
in understanding the mechanism that leads to the steady 
state. This is shown in Fig. |l|. Whereas clearly the input 
current is a random function between and 1 , the output 



current displays sequences of bursts followed by long pe- 
riods of quiescence similar to the one found by Hwa and 
Kardar in the slow driving regime. We also computed 
its power spectrum S(v) (the Fourier transform of the 
output current-current correlation) which appears to be 
white noise in all cases. This is related with the fact that 
our system corresponds to a non-interacting avalanche 
regime in their language fl5|| . 

We now turn to the calculus of critical exponents. 
First we consider the exponent z as defined in (H). This 
was computed by plotting the average duration of the 
avalanches as a function of their characteristic average 
lengths. A binning procedure analogue to the one used 
in JjJ was employed. Plots are shown in Fig. ^. Our best 
fit estimates are 1.34 ± 0.02 and 1.65 ± 0.02 in d = 2, 3 
respectively, compatible with the BTW values which are 
4/3 and 5/3. Remarkably, these results are also in per- 
fect agreement with the RG results of Ref. jlO) which 
are 1.36 (d = 2) and 1.68 (d = 3). The RG analysis 
was performed on a Langevin equation where the driv- 
ing and the relaxation time scales are comparable (and 
hence not well separated). Furthermore the strong (infi- 
nite) non-linearity appearing in the continuum analogue 
of Eq. (|]J) , was regularized and the result was analyzed 
within a one-loop RG scheme. In view of all these ap- 
proximations, the aforementioned closeness in the two 
results is rather surprising. We shall come back to this 
issue later on. 

Another interesting critical exponent is the avalanche 
exponent size r s defined by the relation 



p(S) = S- T °F(S/L*). 



(12) 



Here p(S) is the distribution density of the avalanche 
sizes S, t s is the avalanche exponent and F(x) is a finite 
size function defining the exponent <f> |l(| . The function 
F(x) is assumed to go to a constant for small arguments 
(i.e. large sizes L) and to "regularize" the large avalanche 
behaviour. In order to improve the numerical estimates 
it proves convenient to look at the integrated distribution 
density defined as 



P(S) = / ds p(s). 



(13) 



We have estimated the values of r s in two different ways. 
By plotting the local slope (Fig. 3) and upon a finite 
size procedure (analogue to the one used in Refs. ^ and 
J17J, see Fig. 4). Both procedures yield consistent re- 
sults. In d = 2 our best estimate is 1.288 ± 0.019 which is 
close to the one given in Ref |5) by Liibeck, who reports 
1.282 ± 0.010. It appears however that the two extrapo- 
lations are not identical, since in his analysis the values 
are increasing as L increases rather than decreasing as 
one would expect from a finite size scaling. 

Remarkably, both values are in good agreement with 
the BTW value, thus supporting the claim that the 
Zhang model belongs to the BTW universality class [g]. 
Our d = 3 result is 1.454±0.041 and it supersedes the one 
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reported by Janosi g namely 1.55 which was presumably 
based only on small sizes (and thus too high according 
with our previous discussion). However this disagrees 
with the BTW value 4/3 (see e.g. @) and hence with 
the claim that the BTW and Zhang model belong to the 
same universality class. 

The values of </> were computed from the collapse of 
the curves obtained plotting S Ts P(S) versus S/L^ , that 
is the universal finite size function. We find the best col- 
lapse for 1.80 ± 0.05 and 2.6 ± 0.1. The error bars are 
estimated graphically. A consistent value can be esti- 
mated by plotting the size of the maximum avalanche as 
a function of the size L, which is expected to scale as: 



L*. 



(14) 



A log-log linear fit yields 1.84 ± 0.06 (d = 2) and 2.54 ± 
0.09 (d = 3). A summary of all these critical exponents 
is reported in Tables I and II. 

Let us now turn to the behavior of the roughness as 
defined in Eqs. (|TT|). As mentioned earlier, the dynami- 
cal exponent z can be found from the scaling ansatz (|^) . 
However, as it is usually done in the field of growth pro- 
cesses |l4| , one might think to derive it from the scaling 
of the roughness as well. In Fig || w e plot the rough- 
ness as defined in (|Io|). We find (11a) to hold true with 
fa = 0.282 ± 0.013 and /?/ = 0.391 ± 0.031 in d = 2,3 
respectively. These values were obtained upon using an 
analysis similar to the one exploited to compute r s . From 
the collapse plot one can then infer the value of Zf ap- 
pearing in ( |ll"a| ). We find z f = 1.20 ± 0.05 (d = 2) and 
Zf = 1.4 ± 0.1 (d = 3), which are both lower than the 
corresponding value derived from <^j. Similarly to ( |l4| ) 
we have that 



L Zf 



(15) 



with z f = 1.19±0.04 (d = 2) and z f = 1.34±0.04 (d = 3). 
Commonly the equality z = Zf is tacitly assumed to hold 
and we are not aware of any other examples where this 
point was sufficiently emphasized. A simple argument 
can be given here to explain this discrepancy. In usual 
interface growth phenomena the dynamical exponent is 
measured as the scaling of the saturation time with the 
system length, and this saturation occurs when the cor- 
relation length reaches the system length. In our case, 
both lengths do not scale linearly but as £ ~ L v . Thus 
these exponents need not be identical unless rj = 1. By 
a direct measurement (looking on how the maximum £ 
scales with L) we have found that 77 = 0.922 ± 0.012 and 
i] = 0.897 ± 0.051, for d = 2 and 3, respectively. Accord- 
ing to these scaling arguments we find that the product 
zr\ agrees, within error bars, with the values reported 
for Zf. In certain surface growth models a similar phe- 
nomenon, called anomalous scaling, has been reported 
There it has been observed that the roughening 
exponents are different when measuring the local or the 
global widths. 



Another exponent is derived from the relation Xf — 
(3f Zf which is telling that the roughness, after that the 
avalanche has been completed (i.e. at time T(t)), de- 
creases as L~ Xf . The values Xf — &f z f> according to 
our previous results, are 0.33 and 0.55 in d = 2, 3 respec- 
tively. We now go back to the comparison with the RG 
results. 

As previously hinted, although the exponent z derived 
from (||) is very close to the one derived by RG methods 
on the continuum Langevin analogue of the Zhang model 
||, the (3f and Xf exponents are not. A summary of all 
these values is reported in Table III for compactness. We 
argued previously that this inconsistency is not surpris- 
ing in view of the different physical regimes probed by 
the two cases and of the heavy approximations involved 
in the RG calculation. The apparent equality in the dy- 
namical exponent z then probably hinges on deeper and 
more interesting reasons, and we are planning to consider 
this in a future work. 

Finally, we have also measured the roughness on the 
slow time scale as defined by ([)]). We find that after a 
transient scaling as L d , the roughness tends to a limi t 
which is independent on L (see Fig. ^), i.e. eq. ( p_lb| ) 
holds with /3 S = 0, Xs = and z s = d. 



V. CONCLUSIONS 

In this paper we have studied the infinitely slowly 
driven Zhang model in two and three dimensions. In 
two dimensions this work can be seen as a complement 
of an earlier large sizes study g. On the other hand in 
three dimensions, our results are expected to improve an 
earlier estimate ||. The aim of || was different from 
ours and this could account for the difference. In both 
cases we computed some exponents (notably the </> and 
all the roughness exponents) which were never previously 
considered. Besides being an useful complement to the 
existing literature on the model we also found few un- 
expected results: i) the three dimensional avalanche size 
exponent does not coincide with the BTW value, as the 
two-dimensional value seems to suggest; ii) the exponent 
z computed from the dynamical scaling ansatz does not 
coincide with the one computed from the roughening ex- 
ponent. We have shown that this stems from the non- 
linear scaling of the correlation length £ with the system 
size L; iii) the coincidence between the value of z of the 
Zhang model and the RG value derived on its Langevin 
continuum counterpart, does not extend to other expo- 
nents such as the /3 and the x exponent. 

We believe that all the above issues deserve further 
attention both from the analytical and numerical view 
point. We are currently performing a numerical investi- 
gation on the continuum Langevin equation. This further 
analysis is expected to shed new lights on the approxi- 
mations involved in the RG treatment. 
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FIG. 1. Plot of the total energy J(t), both in (dotted line) 
and out (full line), in d = 2 (a) and d = 3 (b). 
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FIG. 3. Local slope plot for t s as a function of the 
avalanche size S in d = 2 (a) and d = 3 (b). In both cases 
the intermediate most linear part of the largest size was used 
for the computation. 
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FIG. 5. Log- Log plot of Wf(r, L) as function of r for vari- 
ous sizes L. In d = 2 these were L=70 (□), 100 (o), 150 (V), 
200 (A), 300 (o), and in d = 3 they were L=20 (❖), 30 (V), 40 
(A), 50 (o), 60 (□). In both cases the solid line corresponds 
to the value of /?/ reported in Table III. 



FIG. 4. Finite size plot for r s as a function of the inverse 
lattice size 1/L in d = 2 (a) and d — 3 (b). 
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TABLE I. Critical exponents r s and <f) in d — 2, 3. The 
values indicated by (a) and (b) refer to the BTW (l7) and the 
j_respectively. The exponent <f> given here 

!)• 



previous works [||^ 
is computed from (1 



d 




T s (a) 


r.(b) 




P(a) m 


2 


1.288 ±0.019 


1.293 


1.282 ± 0.010 


1.84 ±0.06 


2 


3 


1.454 ± 0.041 


4/3 


1.55 


2.54 ±0.09 


3 



TABLE II. Dynamical critical exponent in d = 2, 3. The 
first column corresponds to (5), whereas the second column 
is computed from (15). Finally the last two columns indi- 
cated by (a) and (b) are the BTW and the RG values Q, 
respectively. 



d 



z f 



z ( a ) 



z(b) 



1.34 ±0.02 
1.65 ±0.02 



1.19 ±0.04 
1.34 ± 0.04 



4/3 
5/3 



1.36 
1.68 



0.10 
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0.00 




TABLE III. Roughness exponents /3 and x m d = 2, 3. The 
values j3f and Xf — Pf z f are computed here while the others 
are the RG values [pi. 



d 


Pf 


P 


X f 


X 


2 


0.282 ±0.013 


-0.26 


0.33 ± 0.03 


-0.36 


3 


0.391 ±0.031 
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0.55 ±0.08 
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FIG. 6. Plot of W 3 (t, L) as function of t for various sizes L 
both for d = 2 (a) and d — 3 (b). The values used for the sizes 
are the same of the previous figure, the larger L the slower 
the growth. 
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